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1  Objectives 

The  overall  goal  of  the  Phase  I  effort  was  to  develop  the  basic  software  framework  for  an 
object-oriented  toolbox  for  modeling,  control  design  and  analysis  of  distributed  parameter 
systems  based  on  PDESolve,  our  C++  class  library  for  simulating  partial  differential 
equation  (PDE)  systems.  Specific  technical  objectives  are  listed  below: 

1.  Develop  the  capability  to  calculate  the  basic  system  matrices  through  spatial  dis¬ 
cretization  of  linear  PDE  control  problems  using  the  finite-element  method. 

2.  Develop  an  interface  with  MATLAB  to  allow  the  user  easy  access  to  the  control 
methodologies  in  MATLAB  such  as  LQG,  Uoo  and  ^-synthesis. 

3.  Develop  the  capability  to  perform  model  order  reduction  at  the  PDE  level  by  using 
an  appropriate  eigenfunction  family  as  the  basis  for  the  spatial  discretization. 

4.  Develop  the  capability  for  simulating  the  open-loop  and  closed-loop  response. 

5.  Demonstrate  the  Phase  I  product  on  a  model  PDE  control  problem. 

2  Work  Performed  and  Results 

2.1  Approach 

The  approach  taken  in  Phase  I  was  to  leverage  PDESolve  to  rapidly  develop  the  ability 
to  perform  PDE-based  control.  PDESolve  is  a  C++  (object-oriented)  class  library  that 
enables  a  high  level  expression  of  ordinary  and  partial  differential  equations,  and  interfaces 
with  real-world  engineering  tools  such  as  CAD  systems.  PDESolve  was  extended  with  new 
objects  that  are  specific  to  control. 

PDESolve  is  a  new  class  of  tool,  offering  mixed  symbolic  and  numeric  computing,  open 
Ch — (-  architecture,  easy  interface  with  engineering  tools  and  the  potential  for  scalable  per¬ 
formance.  Existing  tools  such  as  MATLAB  and  Mathematica  have  some  but  not  all  of  the 
above  characteristics.  The  reusability  and  flexibility  that  come  with  the  object-oriented 
and  math-based  architecture  of  PDESolve  allows  dramatically  faster  development  of  en¬ 
gineering  solutions.  These  solutions  are  often  one  to  two  orders  of  magnitude  shorter  in 
terms  of  line  count  compared  with  programming  in  C  or  FORTRAN,  and  take  accordingly 
less  time  to  develop.  More  than  $4  million  dollars  have  been  invested  in  PDESolve  and 
related  projects  to  date. 

2.2  Results 

In  the  Phase  I  effort,  we  have  developed  the  core  functionality  for  the  control  toolbox  and 
applied  it  to  various  PDE  control  problems,  thus  demonstrating  the  technical  feasibility 
of  the  proposed  Phase  II  effort.  The  main  accomplishments  in  Phase  I  are  listed  below: 

1.  High-level  specification  of  PDE  control  problem 

2.  Calculation  of  system  matrices 

3.  Interface  with  MATLAB  providing  access  to  MATLAB  functionality  from  within 
the  control  toolbox 
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4.  Simulation  of  open  and  closed-loop  response 

5.  Demonstrations  on  model  PDE  problems  in  heat  and  vibration  attenuation 

6.  Model  reduction  at  the  PDE  level  using  the  modal  basis 

Each  of  the  above  items  is  discussed  in  detail  below.  This  discussion  is  illustrated  using  a 
simple  physical  example.  Let  U(t,x,y)  be  the  temperature  distribution  in  a  thin,  homoge¬ 
neous  plate.  The  objective  is  to  regulate  the  temperature  of  the  plate  using  a  distributed 
actuator  u{t,x,y)  that  can  source  and  sink  heat.  The  optimal  control  problem  can  be 
stated  as  follows:  Find  the  control  u  that  minimizes  the  performance  index 

J{u)  =  JQf™0  {W,x>y)\2  +  R\u(t,x,y)\2}dtdxdy  (1) 


subject  to  the  PDE  system 


=  k 


dU 
dt 

y  =  U 


'd2U  d2US 
dx 2  dy2 


+  u 


U(t,0,y)  =  U{t,20,y)  =  U{t,x,  0)  =  U{t,x,  20)  =  0 

.  »  f  2  8  <  x  <  12,  8  <  y  <  12 

U{  ,x,y)  |  o  everywhere  else 

where  y(t,x,y)  is  the  observed  variable  (a  distributed  sensor). 


(2) 

(3) 

(4) 

(5) 


High-level  specification  of  PDE  control  problem 

We  have  developed  the  capability  to  specify  the  PDE  control  problem  at  the  continuous 
level.  Since  we  use  the  finite-element  method  for  spatial  discretization,  the  problem  is 
input  into  the  toolbox  in  the  weak  form  which  for  eqs.  (2)-(3)  is  given  by 

[  ®E.vdO  +  k  [  VU-Vvdtl-  f  ^-vdT-  [  uvdQ  =  0  (6) 

J a  dt  Jn  Jr  dn 

[  {y-U)dO  =  0  (7) 

Jn 

where  v  is  the  test  function.  Since  the  boundary  conditions  are  Dirichlet,  the  underlined 
term  is  evaluated  using  Lagrange  multipliers  and  is  set  to 

f  ^Lvdr=  f  XvdT+  f  SXUdT  (8) 

Jr  dn  Jr  Jr 

SX  and  v  being  linearly  independent  test  functions.  The  above  weak  form  can  be  input 
directly  into  the  toolbox  as  shown  in  the  code  segment  below. 

Function  U(2,  Scalar);  //  Unknown  function 
Function  v(2,  Scalar,  VARIATIONAL);  //  Test  function 
Function  lambda(2,  Scalar);  //  Constraint  dofs 

Function  varLambda(2,  Scalar,  VARIATIONAL);  //  Constraint  test  func 
Function  uControl(2,  Scalar);  //  Control  variable 
Function  y(2,  Scalar);  //  sensor  variable 
Real  kDif fusion  =  0.25; 
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//  Form  cell  complex  and  FE  mesh  using  third-party  mesher . 

CellComplex  cc  =  rectMeshGen(nx,  ny,  0.0,  20.0,  0.0,  20.0); 

FEMesh  mesh(cc,  2); 

//  Weak  formulation  of  the  PDE. 

Equations  w(mesh) ; 

w  =  Integral (v*(dt*U)  +  kDiffusion*(grad*U)*(grad*v)) 

-  Integral (uControl*v) ; 

//Weak  formulation  of  the  sensor  equation. 

Equations  s(mesh); 
s  =  Integral (  v*y  -  v*U) ; 

//  Impose  essential  (Dirichlet)  boundary  conditions  using  Lagrange 
//  multipliers 

w["x=0"]  =  lambda*v  +  varLambda*U; 

w["x=L"]  =  lambda*v  +  varLambda*U ; 

w["y=0"]  =  lambda*v  +  varLambda*U; 

w["y=L"]  =  lambda*v  +  varLambda*U; 

Note  the  one-to-one  correspondence  between  the  weak  form  and  its  specification  in  the 
control  toolbox. 

If  instead  of  the  above  continuously  distributed  control,  we  had  actuators  at  specified 
locations  on  the  plate,  we  would  identify  the  actuator  locations  with  labels  while  gener¬ 
ating  the  geometry /cell-complex.  We  then  apply  the  control  term  only  over  the  region 
represented  by  the  actuator  labels  as  shown  below: 

w  =  Integral (v*(dt*U)  +  kDiffusion*(grad*U)*(grad*v)) 

-  IntegralC'actuator",  uControl*v) ; 

where  the  location  of  heat  source/sink  is  identified  by  the  label  “actuator”. 


Calculation  of  system  matrices 

The  discrete  control  problem  on  the  given  finite-element  mesh  is  formed  using  the 
ControlProblem  class  and  the  system  matrices  are  obtained  using  the  getMatnces 
method  on  it  as  shown  in  the  following  code  excerpt: 

ControlProblem  problem(mesh,  dt,  w,  s,  List (U, lambda) ,  List(v.varLambda) , 

uControl ,  y) ; 

DenseMatrix  A; 

DenseMatrix  B; 

DenseMatrix  C; 

problem. getMatrices (A,  B,  C) ; 

The  arguments  passed  to  the  ControlProblem  constructor  are  the  finite-element  mesh 
mesh,  the  time  derivative  dt  to  indicate  the  time  direction,  the  weak  form  of  the  governing 
equations  w  and  s,  the  list  of  unknown  functions  List (U, lambda),  the  corresponding 
list  of  variational  functions  List(v.varLambda),  the  control  variable  uControl  and  the 
sensor  variable  y.  Most  of  this  is  specific  to  LQR  problems  and  was  done  to  simplify  the 
description  of  the  problem;  a  D  matrix  and  a  partition  of  B  and  C  will  be  necessary  for 
general  control  problems. 
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Interface  with  MATLAB 

We  have  built  an  interface  to  MATLAB  that  allows  easy  access  to  all  of  MATLAB’s 
capabilities  from  within  our  control  toolbox.  Users  can  perform  optimal  control  design 
(LQG,  Roo,  /J-synthesis,  C\  optimal  control  etc.),  plant  and  controller  model  reduction, 
simulate  the  designed  controller  on  the  discretized  plant,  and  in  general,  perform  arbitrary 
MATLAB  operations  based  on  the  system  information  passed  along  through  the  interface. 

The  software  architecture  of  the  interface  is  interprocess  communication  (IPC).  The 
MATLAB  process  and  the  toolbox  operate  as  two  separate  processes,  possibly  running 
on  different  machines.  Data  is  passed  between  them  via  IPC  and  a  C-language  interface 
provided  by  Mathworks  Inc.  Several  C++  wrapper  classes  to  the  C-language  interface 
have  been  implemented,  in  order  to  simplify  use  and  provide  a  consistent  interface  with 
our  data  types.  The  resulting  system  allows  both  applications  to  operate  on  data  in  their 
native  formats,  resulting  in  a  highly  efficient  system. 

In  the  heat  control  example  considered  above,  the  feedback  gains  are  calculated  us¬ 
ing  the  lqry  function  in  MATLAB  which  performs  linear-quadratic  regulator  design  for 
continuous-time  systems.  The  syntax  for  using  the  lqry  function  is 

K  =  lqry (SYS, Q,R) 

where  SYS  is  the  continuous-time  system,  K  is  the  optimal  gain  matrix  such  that  the 
state-feedback  law  u  =  -KU  minimizes  the  cost  function 

J(u)  =  J  {yTQy +  uTRu}  dt  (9) 

In  the  above,  y(t)  and  u(t)  are  spatially  discretized  versions  of  y(t,  x ,  y)  and  u(t,  x,  y).  The 
following  code  segment  illustrates  the  calculation  of  the  gain  matrix  using  the  MATLAB 
interface: 

//Start  a  MATLAB  process  on  remote  machine  and  display  output  on 
//local  machine : 

MatlabShell  matlab("rsh  redfish  matlab5  -display  seneca:0.0") ; 

//Transfer  A  and  B  matrices  to  MATLAB: 

matlab.put(A,  “A”);  //  Pass  matrix  to  MATLAB,  associate  it  with  name  "A" 
matlab.putCB,  “B”);  //  Pass  matrix  to  MATLAB,  associate  it  with  name  "B" 
//Calculate  gain  in  MATLAB. 

mat lab . eval ( " sys  =  ss(A,  B,  C,  0);"):  //Form  continuous-time  system 
int  nSens  =  C.numberOfRowsO ;  //Number  of  sensors 
int  nAct  =  B.numberOfColsO ;  //Number  of  actuators 
char  cmdStr [100]; 

sprintf  (cmdStr ,  "K=-lqry(sys ,  eyeC/.d),  0.  l*eyeC/.d)  )  ; " ,  nSens,  nAct); 
mat lab . eval (cmdStr) ; 
mat lab. eval ("save  K"); 

//Import  gain  from  MATLAB. 

DenseMatrix  gain; 
matlab.  get  (gain,  “K”); 

Note  that  in  matlab . eval ( . . .),  the  part  within  the  parentheses  is  what  the  user  would 
normally  type  in  within  MATLAB.  The  output  weighting  matrix  Q  has  been  set  to  identity 
and  the  control  weighting  matrix  R  has  been  set  to  0.1  times  identity.  Alternately,  an 
interactive  MATLAB  shell  can  be  opened  using 
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matlab . interact () ; 

and  the  feedback  gain  calculated  interactively: 

>  K  =  -lqryCsys ,  eye(nSens),  0. l*eye(nAct)) 

Simulation  of  open  and  closed-loop  response 

We  have  implemented  the  capability  to  simulate  open-loop  and  closed-loop  response  using 
the  openLoopSolve  and  closedLoopSolve  methods  on  the  ControlProblem  class.  The 
use  of  these  methods  for  the  heat  control  example  being  considered  is  shown  below.  Note 
that  the  feedback  gains  are  calculated  in  MATLAB  as  discussed  earlier. 

//Specify  time  stepping: 

TimestepSpec  step(initialTime,  finalTime,  nSteps) ; 

//  Simulate  open-loop  response: 

Array<Function>  solnOpen  =  problem. openLoopSolve (uO,  step, 

BackwardsEuler ( ) ,  Filter (filterFunc)) ; 

//Ouput  in  MATLAB  movie  format: 

output Animate ("openLoop" ,  animateTimeSeries (solnOpen) ) ; 

//Simulate  closed-loop  response  with  the  previously  calculated  gain: 
Array<Function>  solnClosed  =  problem. closedLoopSolve (gain,  uO,  step, 

BackwardsEuler () ,  Filter (filterFunc) ) ; 

//Ouput  in  MATLAB  movie  format: 

output  Animat eO'closedHe at 2D" ,  animateTimeSeries ( solnClosed) ) ; 

//  Output  L2Norm(T)  vs  time. 

ofstream  of ("L2norm.dat") ; 

for  (int  i=0;  i<solnOpen.length() ;  i++)  { 

0f  «  i  «  "  "  «  12Norm(soln0pen[i] )  « 

11  "  «  12Norm(solnClosed[i] )  «  endl; 

> 

The  arguments  to  closedLoopSolve  are  the  gain,  the  discretized  initial  condition  uO, 
the  time-step  specification  step,  the  integration  method  BackwardsEuler  and  the  output 
filter  Filter  (filterFunc).  The  discretized  form  of  the  initial  condition  (5)  is  generated 
with  the  following  code: 

Function  uO (2,  Scalar,  uOFunc) ; 

FEDiscretizer  disc (mesh); 

uO  =  disc. discretize (uO) ;  //Discretize  uO  on  mesh 

with  the  function  definition  for  uOFunc  being  given  by 

Real  uOFunc (const  Coordsfe  x)  { 

if  (  x[o]  >=  8  kk  x[0]  <=  12  kk  x[l]  >=  8  kk  x[l]  <=  12) 

return  2 . ; 
else 

return  0 . ; 

> 
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Natural  BCs  Natural  BCs 


Figure  1:  Schematic  of  the  plane  vibration  demonstration  problem.  The  control  variable 
is  oxy  over  the  diamond-shaped  region  in  the  center. 


The  output  filter  allows  any  general  mathematical  operation  to  be  performed  on  the  so¬ 
lution  field  at  each  time-step  and  only  the  result  from  applying  the  filter  is  stored.  This 
provides  the  user  with  the  flexibility  of  storing  only  pertinent  data  (  say,  for  example,  the 
solution  at  certain  critical  points). 


Demonstration  on  model  problem:  LQR  control  of  plate  vibration 

We  next  demonstrate  our  code  on  the  LQR  control  of  the  vibration  of  a  plate.  We  model 
planar  motion  of  a  thin  square  plate  using  the  plane  stress  approximation.  The  system  is 
controlled  by  a  stress  actuator  on  a  square  at  the  center  as  shown  in  Figure  1;  the  control 
parameter  is  the  xy  component  of  the  stress  at  the  actuator. 

The  weak  form  with  respect  to  a  test  function  v  of  the  governing  equation  is 


J  pv  •  iidfl  +  J  tv  •  u dQ,  -t-  <j  •  SedCl  +  -r) 


d5vx  dvy 
[  dy  dx 


dQ.  =  0 


(10) 


We  impose  the  constraint  u  =  0  at  the  surface  y  =  0,  and  assume  natural  boundary 
conditions  at  all  other  surfaces.  The  plate  is  initially  displaced  in  the  plane  as  shown  in 
Figure  2;  the  initial  velocity  field  is  zero  everywhere. 

Time  series  for  the  open  and  closed  loop  solutions  are  shown  in  Figure  3. 

Fragments  of  the  code  used  to  solve  this  problem  are  presented  below.  For  brevity, 
we  have  edited  out  some  initialization  code,  the  calls  to  the  MATLAB  interface,  and 
PDESolve  postprocessing  steps;  that  code  is  similar  to  that  in  the  heat  equation  example 
above.  Note  that  although  the  dynamic  elasticity  equations  are  considerably  more  complex 
that  the  heat  equation,  the  overall  structure  of  the  code  is  nearly  as  simple  as  in  the  heat 
equation  example. 


//  geometry  initialization  code  deleted  for  brevity. 

//  We  have  created  a  cell  complex  with  labeled  actuators. 
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Figure  2:  Initial  displacement  of  the  elastic  plate.  The  initial  velocity  is  zero. 


time 


Figure  3:  Time  evolution  of  L 2  norm  of  displacement  field  for  open  loop  and  closed  loop 
solutions  of  the  plate  vibration  control  problem. 
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FEMesh  mesh(cc,  2); 


DiffOp  dx(l ,0) ; 
DiffOp  dy(l , 1) ; 
DiffOp  dt(l ,2) ; 
DiffOp  dtt(2,2) ; 


Function  u(2, Tensor (2, 1) ) ; 

Function  v(2,Tensor(2, 1) .VARIATIONAL) ; 
Function  lagr (2, Tensor (2,1)) ; 

Function  var Lagr (2, Tens or (2,1) .VARIATIONAL) ; 
Function  z(2,  Scalar); 


//  displacement  field 
//  test  function 
//  lagrange  multiplier  for  BCs 
//  test  function  for  BCs 
//  control  variable 


Real  damp  =  0.05; 

Real  mu  =  E/(2.0*(l+nu)) ; 

Real  lambda  =  E*nu/((l+nu)*(l-2*nu)) ; 


//  damping  constant 
//  Lame’  constants 


//  form  stress-strain  relation 
Tensor  material (3, 2) ; 
material  =  List(  List(lambda+2*mu, 
List (lambda, 
List(0, 


lambda,  0  ) , 

lambda+2*mu,  0  ) , 

0 ,  mu) ) ; 


Expr  strainDisp  =  List (List (dx, 0) ,  List(O.dy),  List(dy,  dx)); 
Expr  strain  =  strainDisp  *  u; 

Expr  varStrain  =  strainDisp  *  v; 

Expr  stress  =  material  *  strain; 


//  set  up  weak  form  of  equations 
Equations  w(mesh); 

w  =  Integral (v*dtt*u  +  0.05*v*dt*u  +  stress*varStrain)  + 

Integral ("actuator",  v*dtt*u  +  0.05*v*dt*u  +  stress*varStrain  + 

z*(dx*v[l]  +  dy*v[0])); 

//  variational  enforcement  of  dirichlet  BCs 
w["y=0"]=  lagr*v  +  varLagr*u; 


//  assemble  control  problem 

ControlProblem  problem (mesh,  dt,  w,  List(u.lagr) ,  List(v.varLagr) ,  z) ; 


//  form  discrete  initial  conditions 
Function  u0(2.  Tensor (2,1),  initPosFunc) ; 
FEDiscretizer  disc (mesh); 
uO  =  disc .discretize (uO) ; 

Function  vO (2,  Tensor(2,l),  initVelFunc) ; 
vO  =  disc .discretize (vO) ; 
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// 

//  calls  to  matlab  deleted  for  brevity. 

//  see  the  first  example  in  this  section  for  matlab  interface  example  code. 

// 

TimestepSpec  step(0.0,  10.0,  60); 

Array<Function>  solnOpen  =  problem. openLoopSolve (List (uO, vO) , 
step,  BackwardsEulerO); 

Array<Function>  solnClosed  =  problem. closedLoopSolve (gain,  List(uO.vO), 
step,  BackwardsEulerO); 

//  postprocessing  and  animation  calls  deleted  for  brevity 

> 

Model  reduction  at  the  PDE  level  using  the  modal  basis 

The  examples  presented  so  far  have  used  a  finite-element  discretization  of  the  PDE.  That 
is  acceptable  for  control  of  small  demonstration  problems,  but  for  real-world  PDE  control 
problems  the  number  of  degrees  of  freedom  becomes  prohibitively  large.  A  useful  alter¬ 
native  is  to  represent  the  system  in  terms  of  a  small,  dynamically  meaningful  set  of  basis 
functions  such  as  an  eigenmode  decomposition  or  a  proper  orthogonal  decomposition.  De¬ 
termination  of  an  appropriate  set  of  functions  through  an  eigenmode  analysis  or  POD  is 
itself  a  computationally  intensive  problem,  but  once  done,  it  allows  the  control  designer 
to  work  with  a  far  simpler  system. 

We  have  implemented  the  infrastructure  needed  to  perform  dimension  reduction  via  a 
modal  representation.  In  the  example  below,  we  consider  control  of  the  heat  equation  on  a 
2D  square.  Control  is  accomplished  by  four  point  source  actuators  inside  the  square.  We 
represent  the  solution  by  the  16  lowest  frequency  modes  on  that  domain;  we  obtain  those 
modes  though  a  finite-element  analysis  using  PDESolve.  Using  the  same  finite-element 
mesh  to  do  the  controls  problems  results  in  a  system  with  225  degrees  of  freedom;  the 

modal  representation  has  reduced  that  to  16. 

The  changes  in  user  code  needed  to  switch  from  a  full  FE  model  to  a  modal  represen¬ 
tation  are  minimal.  As  shown  in  the  code  fragment  below,  essentially  the  only  changes  are 
(a)  to  use  a  SpectralMesh  rather  than  an  FEMesh  object,  and  (b)  the  boundary  conditions 
are  built  into  the  modes  and  are  thus  not  needed  in  the  modal  calculation.  The  boundary 
conditions  were  of  course  used  in  the  computation  of  the  modes. 


//  initialization  omitted;  mode  computation  shown  below 


SpectralMesh  mesh (modes) ; 

Equations  w(mesh); 

w  =  Integral (v* (dt*u)  +  (grad*u)*(grad*v))  +  IntegralC'actuator" ,  z*v)  ; 
ControlProblem  problem(mesh,  dt,  w,  u,  v,  z) ; 

//  matlab  calls,  ODE  solution,  and  postprocessing  omitted 
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Figure  4:  Comparison  of  full  FE  and  reduced-dimension  models  of  heat  equation.  Plotted 
are  times  series  for  L2  norm  of  temperature  for  open  loop  and  closed  loop  as  solved  with 
FE  and  modal  methods. 

In  Figure  4  we  show  results  for  both  the  reduced-dimension  model  and  the  full  FE 
simulation  of  the  point-actuator  controlled  heat  equation  example.  The  open-loop  solu¬ 
tions  are  indistinguishable  in  the  plot.  The  closed-loop  solutions  with  the  two  methods 
compare  well;  the  differences  are  explained  by  the  inability  of  the  truncated  modal  basis 
to  exactly  represent  the  effect  of  point  source  terms. 

We  next  present  a  code  sample  showing  the  mode  computation  in  PDESolve.  The 
PDESolve  code  sets  up  a  general  linear  eigensystem  that  is  then  solved  using  the 
ARPACK  family  of  iterative  eigensolver  routines.  This,  together  with  the  control  code, 
totals  to  under  one  hundred  lines  of  user  code  for  specification  and  solution  of  a  dimension- 
reduced  PDE  control  problem. 

//  initialization  code  omitted 

DiffOp  dx(l,0) ; 

DiffOp  dy(l,l) ; 

Expr  grad  =  List(dx.dy); 

Function  u(dim, Scalar) ; 

Function  v(dim, Scalar .VARIATIONAL) ; 

Function  invAlpha(dim,  Scalar); 

Function  lambda(dim,  Scalar) ; 

Function  varLambda (dim, Scalar .VARIATIONAL) ; 

FEMesh  feMesh(cc .order) ; 

Equations  heatModeEqn(feMesh) ; 
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heatModeEqn  =  Integral (invAlpha* (grad*u) * (grad*v)  +  u*v)  ; 
heatModeEqn["x=0"]  =  varLambda*u  +  lambda*v; 

heatModeEqn ["x=L"]  =  varLambda*u  +  lambda*v ; 

heatModeEqn["y=0"]  =  varLambda*u  +  lambda*v; 

heatModeEqn ["y=L"]  =  varLambda*u  +  lambda*v; 

EigenProblem  eigen (mesh,  heatModeEqn,  List(u,  lambda), 

List(v,  varLambda) ,  invAlpha,  ConstrainedARPACKSym(EIG_LARGEST,  16)); 

Array<Function>  modes (0); 

Array<Real>  eigenvalues (0) ; 

eigen. solve (eigenvalues,  modes); 

3  Personnel  Involved 

Dr.  Gal  Berkooz 
Dr.  Rajesh  Bhaskaran 
Dr.  Kevin  Long 

4  Publications 


None. 


